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Abstract: We introduce a new computational approach for femtosecond 
pulse propagation in the transparency region of gases that permits full 
resolution in three space dimensions plus time while fully incorporating 
quantum coherent effects such as high-harmonic generation and strong-field 
ionization in a holistic fashion. This is achieved by utilizing a one- 
dimensional model atom with a delta-function potential which allows for a 
closed-form solution for the nonlinear optical response due to ground-state 
to continuum transitions. It side-steps evaluation of the wave function, 
and offers more than one hundred-fold reduction in computation time 
in comparison to direct solution of the atomic Schrodinger equation. To 
illustrate the capability of our new computational approach, we apply it to 
the example of near-threshold harmonic generation in Xenon, and we also 
present a qualitative comparison between our model and results from an 
in-house experiment on extreme ultraviolet generation in a femtosecond 
enhancement cavity. 
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1. Introduction 

The goal of this paper is to present a new computational approach for nonlinear light-matter 
interactions in atomic and molecular gases for incident femtosecond pulses whose spectra lie 
in the transparency region, the approach opening the door to fully resolved propagation in 
three spatial dimensions plus time, or (3+1) dimensions. For sufficiently low input intensities 
when ionization of the constituent atoms is negligible and the atoms remain dominantly in their 
ground state, suitable space-time resolved computational approaches already exist and incorpo- 
rate the linear dispersive properties of the medium, both refractive-index and absorption, along 
with a model for the bound state electronic nonlinearity, typically a Kerr-type nonlinearity 
possibly with saturation, and also a time-delayed nonlinear response to model Raman effects 
for molecules 0"). For higher intensities when atomic ionization becomes relevant a number 
of effects turn on that involve ground state to continuum state transitions. These transitions 
can greatly modify field propagation via nonlinear generation of freed electrons and associated 
nonlinear losses and refraction, and high harmonic generation (HHG). Typically the nonlinear 
optical response arising from ground state to continuum transitions are dealt with in a piecemeal 
fashion. For example, the polarization source term for high harmonics is often calculated using 
the strong-field approximation [2| whereas the rate of ionization is treated in a separate calcu- 
lation using Keldysh theory or one of the subsequent variants such as ADK theory (see for a 
useful compendium of different ionization-rate formulas). Such an approach has met with con- 
siderable success to date for simulations of HHG (ref. J4| has a readable yet detailed description 
of the method) and also light filament propagation in gases (TJ, but it would clearly be advanta- 
geous for the advancement of these fields if the HHG and all ionization related effects could be 
modeled in a holistic manner so that as parameters are varied the relative roles of high-harmonic 
generation and ionization related effects are microscopically consistent without having to tune 
the parameters of separate models. An obvious candidate for such a quantum model is to solve 
the three-dimensional (3D) atomic Schrodinger equation |3]|6) to obtain the quantum averaged 
dipole moment at each point in space, and use this as a source in Maxwell's equation for the 
field propagation. However, due to the restrictive computing requirements needed to solve the 
3D Schrodinger equation it seems clear that this will not be a viable option for the foreseeable 
future. 

The quantum model we employ in the current work is one-dimensional, and it is perhaps 
the simplest system which has both the energetic continuum and a bound state — the minimal 
set of ingredients crucial to describe light-matter interactions at femto-second time scales and 
ionization-level intensities. However, the model has the virtue that it can be implemented to ef- 
ficiently simulate space-time resolved field propagation in (3+1) dimensions. To the best of our 
knowledge Christov J7) was the first to perform simulations of femtosecond pulse propagation 
in one space dimension plus time, or (1+1) dimensions, that incorporated the medium response 
via direct numerical solution of a time-dependent ID atomic model. In a series of recent papers 
Bandrauk and co-workers J8| have presented simulations of ID quantum models coupled to 
field propagation in both (1+1) and (3+1) dimensions, and addressed a range of topics of cur- 
rent interest including HHG and filamentation in gases. One-dimensional atomic models have 
played an important role as model systems for elucidating the physics underlying nonlinear 
light-matter interactions. For example, Su and Eberly J5] introduced the quasi-Coulomb po- 
tential in ID as a model for multi-photon physics. Here we use the 5-potential atomic model. 
This model has one free parameter, the ionization potential, and its spectrum consists of one 
bound state, the ground state, plus the continuum, so it is an ideal model system in which to ex- 



plore the role of ground state to continuum transitions. Indeed, this system has been studied and 
applied in various contexts for many years. The attraction of the model as a toy system to inves- 
tigate various aspects of dynamics in quantum systems stems from the fact that many quantities 
can be obtained exactly 11 1 Oil 111 . Its generalizations to higher dimensions lfl2l . for more com- 
plex potentials 11311141 . and other dynamic equations [15] also exist. Besides its utility in the 
AMO and strong field physics |[T6ljT9l , the ID 5-potential has also been applied in semicon- 
ductors 11 10131 . and used as a test system for numerical and approximate methods [20- 22) . We 
have applied the 5 -potential model to explore the relative roles of the higher-order Kerr effect 
and plasma induced defocusing in filament propagation (23.. 24] . Furthermore, it has recently 
been shown by two of the present authors (MK and JMB) that the 5 -potential model yields 
a closed form solution for the quantum averaged nonlinear current for this model atom in an 
arbitrary time-dependent light field |25|. (The code computing the nonlinear current is avail- 
able as open-source upon request to the corresponding author (MK).) This exact solution side 
steps the need for direct evaluation of the time-dependent atomic wave function and therefore 
represents an enormous savings in computational effort and time. For example, one does not 
need to worry about griding issues for the solution of the ID wave equation and the spurious 
effects of numerical boundary conditions. Our new approach is computationally efficient since 
it employs the closed form for the nonlinear current as a source in the (3+1) field propagation 
equation, and we find gains of more than one hundred in computation time in comparison to 
the approach using integration of the Schrodinger equation. 

The remainder of this paper is organized as follows: In Sec. 2 we describe our computa- 
tional approach and discuss the model equations. In Sec. 3 we include numerical simulations, 
and for this initial paper we utilize near-threshold HHG in Xenon as an illustrative example. 
We chose near-threshold HHG, with concomitant lower-orders in the harmonic spectrum, to 
showcase the capability of the model to cope with an extreme spectral bandwidth in a com- 
pletely unified treatment. Unlike the strong-field approximation, our approach does not rely 
upon identification of the quantum paths that dominate the HHG spectrum, and is valid also for 
very low frequencies. The first part of Sec. 3 describes simulations of HHG in a Xenon gas jet 
to demonstrate that our model can deal with different incident focusing conditions and hence 
phase-matching regimes, and that it reproduces expected qualitative features. In this respect 
we do not claim these results are new, particularly since corresponding simulations of HHG in 
(3+1) already exist that employ the strong-field approximation, but rather are intended to high- 
light the capabilities of the approach. Sec. 3 also includes what we believe are new results to 
model near-threshold HHG in a femtosecond enhancement cavity (fsEC). Here a new feature is 
that the HHG occurs in the presence of the residual plasma that accumulates within the Xenon 
gas jet due to the circulating pulse, and we compare the qualitative features of the harmonic 
spectrum with that of the in-house experiment that has recently been used to produce record 
power levels of extreme ultraviolet (XUV) ll26l . 

2. Computational approach and model 

To contrast the model of this paper to the current state of the art in computer simulation of high- 
harmonic generation and/or femtosecond filamentation, we briefly recall the main features of 
the standard approach. The propagation of the optical and high-harmonic fields are described 
in one-way (i.e. directional) propagation equations, into which the various medium responses 
are coupled. For this purpose, we use the Unidirectional Pulse Propagation Equation (UPPE) 
solver. In the spectral domain the UPPE for propagation along the z-axis takes the general 
form (23 

^ ± (©, Z )=^%(fi), Z ) + ^^P fc± (fi), Zj {E})- 5 -^4 i (fi), Z) {£}), (1) 



where the z-component of the optical wavevector is given by 



k z ((0,k ± ) = ^fi>2(l+ Z (a>))/c2-*2. (2) 

The UPPE describes the evolution of the spatially resolved spectral amplitudes of the scalar 
electric field (fi),z), where co is the frequency, and k± denotes transverse wavenumbers, 
so that the space-time dependent electric field E(r,z,t) may always be reconstructed using a 
3D Fourier (or Hankel) transform over k±,(0. Note that the assumption of a linearly polar- 
ized field is consistent with our adoption of a ID atomic model. Nonlinear effects enter the 
propagation equations through either the polarization (P) or current-density (J) terms. These 
are both evaluated as functionals of the electric field history E(t) at each spatial point, and then 
Fourier-transformed to the spectral representation, e.g. Jk ± (co,z, {E}), that appears in the above 
propagation equation. The linear wavelength-dependent properties of a medium are represented 
through the susceptibilityj(a)). 

In the standard model, there are a number of independent contributions both in the polariza- 
tion and in the current density. They include the optical Kerr effect, ionization in strong fields 
and a freed-electron density evolution equation, separately modeled losses due to ionization, 
avalanche ionization, de-focusing effects of freed electrons, and losses due to freed electrons 
described in terms of an effective Drude-plasma model (TJ. 

Our computational approach for simulating femtosecond nonlinear light-matter interactions 
reflects the fact that the above mentioned effects are all manifestation of the single electronic 
system response to the strong optical field. The goal is to achieve a unified description, and re- 
duce multiple independent model parameters. In general, our model involves three components 
for describing an atomic gas: 

• A description of the linear dispersive properties of the gas via a complex-valued fre- 
quency dependent susceptibility %((0). 

• The ID 5 -potential quantum model for modeling the nonlinear current (J) due to ground 
state to continuum transitions. This will incorporate the effects of ionization, HHG, and 
ionization induced absorption and refraction in a holistic fashion. 

• A Kerr-like nonlinearity and associated nonlinear polarization (P) to capture the nonlin- 
ear optical response due to the ground state to excited bound state transitions. This will 
contribute processes such as four-wave mixing, self -phase modulation, and self-focusing. 
(For a molecular gas we may also add a time-delayed nonlinear response to capture Ra- 
man effects) 

For this initial exposition we shall give a description of each of the three components above 
in the subsections below as appropriate to our illustrative example of near-threshold HHG in 
Xenon gas. The three medium components are coupled into the field-propagation equa- 

tion which encompasses all frequency components from the fundamental to high harmonics. 
While we use the Unidirectional Pulse Propagation Equation, we emphasize that this approach 
may be implemented with any pulse-propagation simulator that resolves the carrier wave of the 
optical field. 

Linear dispersive properties 

Inclusion of the linear dispersive properties of the medium over the full span of harmonics is 
important to properly incorporate the phase-matching and re-absorption that affect the spatial 
evolution of the generated harmonics. In our spectral solver this does not present a problem 
as long as suitable data is available for the index of refraction and absorption properties in a 
sufficiently wide spectral range. The medium is characterized in terms of the linear complex 



susceptibility #(©), and the pulse propagation method utilizes this information in the prop- 
agation constant k z ((0,k±) in (f2| at each frequency or wavelength resolved by the numerical 
simulation. 
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Fig. 1 . Real (dark lines) and imaginary (red line) parts of the linear susceptibility of Xenon 
as employed in the numerical simulations in this paper. The vertical dashed lines mark the 
fundamental wavelength of the fundamental field at 800 nm, and several of its harmonic 
are also indicated. 



Throughout this paper we use Xenon as a representative atomic medium to showcase the 
modeling capabilities, and Fig. Q] shows the real (solid dark line) and imaginary (red line) 
parts of the susceptibility as functions of the angular frequency of the electromagnetic field. 
To create the data set displayed in Fig. Q] we have combined Xenon data downloaded from 
http://henke.lbl.gov/opticalxonstants/ with the refractive-index parametrization obtained from 
Ref. l28ll to create a tabulated representation of the linear susceptibility over a range of fre- 
quencies that spans the fundamental at 800 nm wavelength up to its 35 th harmonic. Note that 
with the real and imaginary parts of the susceptibility function coming from separate sources, 
and with no absorption data available for photon energies below 10 eV, our parametrization 
does not satisfy the Kramers-Kronig relations. To our best knowledge a more complete data 
set is currently not available. Nevertheless, while not a perfect model of the Xenon gas for all 
frequencies this data set employed serves a an illustration example in which all electromagnetic 
frequencies are treated on the same footing by the propagation solver, which operates on the 
full electric field. Generally the accuracy of the linear propagation will be limited by the quality 
of the available data for the linear susceptibility. 

Quantum atomic model for the nonlinear response 

As alluded to in the introduction our computational approach employs a ID quantum model 
with a short-range 5 -potential for the nonlinear current due to ground state to continuum tran- 
sitions. More specifically the corresponding time-dependent Schrodinger equation in atomic 
units takes the form 



l -d}+B8{s)-sF{x) 



iK*,T)=0. (3) 



Here s and x are the space and time variables in atomic units, the strength of the attractive 
potential is governed by the parameter B, with the ionization potential of the ground state equal 
to B 2 /2, and F(t) denotes the time-dependent external field in atomic units applied to the atom. 



The above Schrodinger equation would have to be solved numerically at each spatial point 
resolved by the optical pulse simulator. Such an approach, although feasible for this specific 
quantum system, is very demanding on computing resources. Fortunately, we do not need the 
wave function because it is only the time-dependent dipole moment or current that enters as an 
input in the UPPE (fTJ, and we have recently shown how an exact, closed form solution can be 
obtained for both for an arbitrary external field l25l . More specifically, the solution provides a 
means of evaluating the expectation value for the current in atomic units 



without having to go through solving Eq. [3] for the wave function l//. Rather, the current / 
required for the UPPE is calculated directly, completely by-passing the wave function, which 
amounts to a tremendous savings in computing time while retaining the full coherent quantum 
dynamics of the atomic system as is pertinent to the nonlinear field propagation. Moreover, 
it is possible to exactly eliminate the current component which is linear in the driving field 
strength ll25l . The quantum model can therefore be restricted to contribute only to the nonlinear 
response, which is advantageous when it is combined with the linear dispersive medium. For 
the sake of completeness the formulas for the nonlinear quantum current are summarized in 
the Appendix. The reader is referred to 11251 for details concerning the practical numerical 
implementation. For those interested in using this model in their simulations, the corresponding 
software is available as an open-source upon request to the corresponding author (MK). 

The integration with the pulse propagator solver is in principle the same as for any other 
nonlinear medium response. Having calculated the evolution of the optical field up to a given 
propagation distance, the history of the electric field at a given point in space is converted to 
the external field F(t) in atomic units. This drives the quantum system, and the induced current 
is computed as outlined in the Appendix. The nonlinear current is next converted from the 
atomic units, and multiplied by the number density atoms at the point in space. The resulting 
macroscopic current density is included in the right-hand-side of the UPPE equation. We remark 
that since we only incorporate the nonlinear current from the quantum model into the UPPE, 
we do not double-count by erroneously including linear properties from the ID atomic model. 
Also note that it is not an option to retain the linear part of the model's response instead of 
%{(0) introduced in the previous subsection. This is because in a real medium %((0) originates 
in virtual transitions among a large number of states, and it would be difficult to model this from 
first principles. The linear susceptibility arising from our ID atomic model is far too simplistic 
to capture the chromatic properties of a gas to any realistic degree. 

Next we would like to illustrate that our ID 5-potential atomic model displays features for 
HHG that are qualitatively similar to those obtained using the strong field approximation ap- 
plied to the more exact 3D Hydrogen-like atomic model. To this end Fig.[2]shows an example 
of a harmonic spectrum of the nonlinear current induced in the ID model atom for an ioniza- 
tion potential of 12 eV, characteristic of Xenon, and a ten-cycle pulse of center wavelength 
A = 800nm and peak intensity 1.5 x 10 18 W/m 2 . This harmonic spectrum exhibits the charac- 
teristic high-harmonic generation plateau with a high frequency cut-off. The cut-off predicted 
by the formula from the standard HHG model [f2][29][30) is indicated by the arrow in Fig. [2] for 
our parameters. Furthermore, the harmonic spectrum includes the 9 th harmonic which occurs 
near the ionization potential up to around the 40 th harmonic, and thus constitutes an example 
of near-threshold HHG where the generated harmonics straddle the ionization energy. What is 
shown in Fig. [2] is the spectrum of the nonlinear polarization P that appears in the Maxwell 
equations and not the spectrum of the radiation actually generated. Importantly the harmonic 
spectrum shown is exact, within the context of our ID model, over the whole frequency range 
and in particular at low frequencies. 
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Fig. 2. High-harmonics in the spectrum of the dipole moment induced by a pulse at 
A = 800nm. The intensity was 1.5 x 10 18 W/m 2 , kept constant over duration of ten opti- 
cal cycles. The arrow marks the location of cut-off energy calculated for these conditions. 



Finally we point out some pros and cons of our ID model versus the strong-field approxima- 
tion. The strong field approximation is often utilized in numerical atomic simulations of HHG, 
and describes the 3D atomic system in terms of a single bound state, the ground state, and a 
continuum of free electron continuum states [2]. In this sense it addresses a similar spectrum of 
electronic states as our ID model, albeit in 3D. A distinct advantage of the strong field approx- 
imation is that it can in principle be used for an elliptically polarized driving field, whereas our 
ID model is restricted to linear polarization. On the other hand, ours is an exact solution of a 
well-defined system, and as such it is valid throughout the whole frequency spectrum. Unlike 
the strong field approximation, it accounts for "all electron trajectories" not only those that give 
rise to the harmonic radiation. In particular, the low-frequency components of the nonlinear 
current response affect the propagation of the driver pulse through ionization, de-focusing by 
freed electrons, ionization losses, and the nonlinear focusing. 

The nonlinear Kerr effect 

The 5 -potential atomic model only incorporates the nonlinear contribution of ground state 
to continuum transitions. To capture the nonlinear contribution of virtual ground state to bound 
state electronic transitions we include a term representing the nonlinear Kerr effect of the form 

P{t)=2e Q n b n 2 E 2 (t)E(t), (5) 

where «2 is the nonlinear index and rib is the linear refractive index. This effect enters the 
propagation equation via the nonlinear polarization term, and gives rise to self-focusing and a 
cascade of lower-harmonic radiation. 

3. Illustrations involving near-threshold high-harmonic generation 

HHG in a Xenon gas jet 

The first illustrative example involves near-threshold HHG in a Xenon gas jet. The incident 
pulse in each case is Gaussian in space and time with pulse duration t p = 85 fs and center 
wavelength of 800 nm, and focused spot size wq. We model the Xenon gas jet as a vertical 
column which is tapered along the propagation or z-axis according to the jet pressure profile 
indicated by the dotted line in Fig. [3] This pressure or density profile has a constant region of 



length L = 200 /im between z = 100 — 300 |im and tapers off on each side of this region. For the 
examples in this paper we used 20 Torr for the maximal pressure in the gas jet. Furthermore we 
consider the two distinct focusing cases where the input beam is focused at z = 100 /im close to 
the entrance to the gas jet, and also the case that the input beam is focused at z = 300 jim close 
to the exit of the gas jet (We use this descriptive nomenclature of exit and entrance to the gas 
jet with the caveat that HHG can arise before and after the exit and entrance due to the tapers 
in the gas density). The variation of the on-axis intensity along the z-axis is illustrated in Fig. [3] 
for the cases of a) a narrow beam with wq = 7 jj.m, and b) a wide beam with wq = 15 /im. 
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Fig. 3. For the sake of comparison, we examine four cases of focusing geometry. Here, 
A and B denote geometric focus after and before the center of the gas jet. The nominal 
maximal intensity is kept constant. 



In what follows we shall present examples of fully resolved simulations of pulse propaga- 
tion and harmonic generation. As a means to present our results we have chosen to plot the 
angularly resolved spectra for selected harmonics as they propagate through the gas jet. Such 
spectra have previously been utilized in both numerical simulations and were measured experi- 
mentally 131II32I . For example, with reference to Figs. 4,6,7 each panel shows a color encoded 
plot of the spectral intensity (red being maximum and blue minimum) as a function of the 
transverse wavenumber representing angular spread along the horizontal axis, the center cor- 
responding to zero propagation angle, and frequency along the vertical axis for the harmonic 
indicated. In particular, the frequency axis in each panel is centered on the harmonic order and 
has a spread of one half of the harmonic spacing. Each set of results is arranged into a (3 x 3) 
array of panels, the panels are arranged vertically according to the propagation distance into 
the gas jet, the top panel being the entrance and the bottom panel the exit of the gas jet, and the 
panels are arranged horizontally from left to right for harmonic orders 9,1 1 and 13. 



Tighter focus case 

For the first example we chose a focused spot size of wq = 7 jj.m giving a Rayleigh range of 
zo = ttwq/X = 192 /im for the input field. In this case L ~ zo and we expect that the focusing 
of the beam will impact phase-matching of the HHG via the Gouy phase-shift incurred by the 
Gaussian beam. The calculated angularly resolved spectra are shown in Fig.|4]for (a) focusing 
at the exit, and (b) focusing at the entrance of the gas jet, and we note that there is a marked 
difference between the angularly resolved spectra for the two cases. This is consistent with 
general expectations based on phase-matching |33- 35). For example, on-axis phase-matching 
requires cancellation between the accumulated Gouy phase-shift and the accumulated atomic 
phase-shift due to the quantum path taken by the electron. The Gouy phase-shift is a positive 
quantity whereas the quantum phase is proportional to the z-derivative of the incident beam 
intensity due to the dependence of the quantum phase on the electron quiver energy, which is 
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Fig. 4. Angularly resolved spectra of the 9th, 11th, and 13th harmonics at three differ- 
ent propagation distances through the gas jet for Wq = 7 flm. Each panel shows a fre- 
quency region (vertical axis) corresponding to one half of the harmonic order, centered 
at the given harmonic frequency. The horizontal extent of each panel corresponds to the 
transverse wave-number and thus represents the angle of propagation. 
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b) beam focus before the gas-jet center 
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Fig. 5. Evolution of the conversion efficiency versus propagation distance in the gas jet for 
two different beam focusing conditions. The gas density is a smooth flat-top profile with 
a roughly constant pressure region between 100 and 300 micron. The pulse intensities are 
comparable in both cases, with their maxima approximately located at z = 300 and z = 100 
micron in case a) and b), respectively. 



positive for case (a) with the beam focused at the exit, and negative for the case (b) with the 
beam focused at the entrance to the gas jet. Thus on-axis phase matching is not a possibility for 
the case (a) with the input beam focused at the exit, whereas phase-matching is a possibility in 
case (b) with focusing at the entrance to the gas jet. The difference in the angularly resolved 
spectra evident for cases (a) and (b) in Fig. |4]may thus be traced to expected differences in 
phase-matching conditions 041 . 

If we further examine the angularly resolved spectra at the exit of the gas jet, shown as 
the lower set of panels in Fig.|4]for cases (a) and (b) we observe further consequences of the 



different phase-matching conditions: For case (a) for which strict on-axis phase-matching is not 
possible we see that the off-axis HHG emission is more significant with respect to the on-axis 
emissions for each harmonic order than for case (b) where on-axis phase-matching is possible. 
This is particularly pronounced for the 13 harmonic where case (a) is dominated by off-axis 
emission whereas case (b) is dominantly on-axis emission. So the results of our simulations are 
compatible with expectations based on on-axis phase-matching 0311341 . 
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Fig. 6. On-axis and off-axis propagating harmonic components exhibit different phase- 
matching behavior. The top panel show the angularly resolved spectra of the 9th harmonic 
at the beginning, in the middle and at the end of the gas jet corresponding to case (b) 
from Fig. [4] Energy accumulates preferentially in the off-axis (conical) components during 
the propagation in the center of the jet. Later, the on-axis component becomes relatively 
stronger, as the pulse starts to exit from the jet and the gas pressure it experiences starts to 
decrease. This behavior correlates with the evolution of the energy accumulated in the 9th 
harmonic (shown as the red line in the lower panel). Note that this is the harmonic which 
suffers the strongest losses. 

To continue our discussion Fig.[5]shows the evolution of the conversion efficiency for a vari- 
ety of marked harmonics during the propagation through the gas jet. (The conversion efficiency 
is obtained by integrating the angularly resolved spectrum over the transverse wavenumber 
plane for a given harmonic order.) It compares the cases of (a) focus at the exit (left), and focus 
at the entrance to the gas jet (right), and shows that the efficiency is higher in the latter case, 
as is generally accepted to be the case ll4l|34l. Furthermore, Fig.|4]reveals that the evolution of 
the harmonics during the pulse propagation through the jet is distinctly non-monotonic, that is, 
different angular portions of the harmonic fields dominate at different propagation distances. 
To illustrate that our simulation capability can capture the complex interplay between HHG 
processes and phase-matching considerations due to variation in the incident field and gas jet 
pressure, Fig.|6]plots the evolution of the harmonic conversion efficiency for several harmonics, 
the upper panels showing the evolution of the angularly resolved spectrum for the 9 th harmonic 
at three propagation distances. At the shortest and longest propagation distances note that the 
spectrum has its maximum on-axis, whereas at the on-axis and off-axis become comparable 
in the middle. It is interesting that the turnaround in the conversion efficiency for the 9th har- 
monic that occurs for a propagation distance of around 200 jJ.m is therefore dominated by a 
reduction in the off-axis as opposed to the on-axis emission. This highlights the point that the 



spatial evolution of the various harmonics involves a complex interplay between on-axis and 
off-axis emissions, and it is not generally sufficient to assess the conversion efficiency based on 
the on-axis phase-matching considerations. 

Weaker focus case 

For this case we chose a focused spot size of wo = 15 fim giving a Rayleigh range of zo = 
tiwq/X = 883 jxm for the input field. Now, L < < zo and we expect that the focusing of the beam 
will not impact phase-matching of the HHG in a major way. The calculated angularly resolved 
spectra for this case are shown in Fig.|7]for the cases of (a) focusing at the exit, and (b) focusing 
at the entrance of the gas jet, and we note that there is no major differences difference between 
the angularly resolved spectra for the two cases. There are small quantitative differences though, 
in particular in the visibility of the ring structure (cf. the 1 1th harmonic far-field patterns). 



a) beam focus after gas-jet center b) beam focus before the gas-jet center 




Fig. 7. As expected, for a wider beam (wo = 15/lm), the precise location of the beam 
focus is less important. Angularly resolved spectra are shown for the 9th, 11th, and 13th 
harmonics at three different propagation distances through the gas jet. Each panel shows a 
frequency region (vertical axis) corresponding to one half of the harmonic order, centered 
at the given harmonic frequency. The horizontal extent of each panel corresponds to the 
transverse wave-number and thus represents the angle of propagation. 

Figure [8] provides a comparison of the conversion efficiency evolution versus the propaga- 
tion distance for the two focusing cases. This figure shows that while the optimal position of 
the beam focus is still at the entrance into the gas jet, as generally accepted 1341 . for weaker 
focusing of the input field the qualitative behavior is relatively insensitive to the position of the 
input beam focus with respect to the gas jet. 

Finally we would like to point to one other feature from our simulations that correlates with 
previous studies: Inspection of the high resolution image for the angularly resolved spectra for 
the 1 1th harmonic in Fig.Qreveals that the output spectra in the bottom row display ring struc- 
tures. Such structures have previously been seen experimentally and appeared in propagation 
simulation in (3+1) dimensions employing the strong-field approximation 03111341 . Moreover 
the appearance of these angle-frequency spectral features has been associated with the inter- 
ference between the short and long path contributions to HHG. Recently, Teichman et al. fl32l 
demonstrated, experimentally and numerically, that rings in the angle-frequency spectra can be 
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Fig. 8. Evolution of the conversion efficiency versus propagation distance in the gas jet for 
two different beam focusing conditions. The gas density is a smooth flat-top profile with 
a roughly constant pressure region between 100 and 300 micron. The pulse intensities are 
comparable in both cases, with their maxima approximately located at z = 300 and z = 100 
micron in case a) and b), respectively. 



attributed to the interferences in the single-atom response in the transverse cross-section pro- 
file, and that their appearance depends on particular conditions such as intensity and specific 
harmonics. It is satisfying that our computational approach can reveal such structures without 
resort to the strong field approximation. Moreover, inspection of near-field cross-sections (not 
shown) of the frequency-filtered atomic response shows that complex far-field patterns correlate 
with the complex spatial transverse distribution of the HHG "source." This suggests, in keeping 
with ll32l . that dynamics which controls the transverse properties of the atomic response play 
equally important role in macroscopic selection through phase matching. 

HHG in a femtosecond enhancement cavity 

As a final illustration, we show results for harmonic generation in a xenon gas jet placed 
inside a passive femtosecond enhancement cavity (fsEC). This experimental configuration has 
been demonstrated as a means to generate XUV frequency combs with the harmonic light being 
generated at each cavity round trip 13611371 . In this geometry the ionized gas jet behaves as a 
nonlinear medium at the focus of a 6 m ring cavity. A pair of 15 cm radius of curvature focusing 
mirrors lead to an intracavity spot size of wq = 15 fim within the jet. The fsEC has a 1% input 
coupler that allows pulse energy enhancement over 200 times relative to the incident pulse train 
leading to peak intensities in excess of 1 x 10 14 W/cm 2 when the gas jet is off. The harmonics 
are coupled out of the cavity using a thin sapphire plate, aligned at Brewster's angle for the 
fundamental pulse, and resolved spatially by reflection from an XUV grating. 

For initial conditions in these simulations, we have utilized a previously calculated temporal 
profile of the pulse just before it enters the gas jet. This pulse profile was obtained from a ID 
steady-state calculation of the fundamental pulse building up inside the fsEC in the presence of 
the xenon gas and includes the effects of cavity dispersion. The high nonlinearity of ionization 
leads to chirped pulses circulating in the cavity and limits the achievable intensities in this 
geometry. Details of this model and its implications for HHG in a fsEC can be found in Ref. 
J38). 

Yet another modification of the current model consists in including the surviving plasma 
in the jet based on the estimated levels calculated in i38l . Because the 20 ns cavity round- 
trip time is too short for the plasma to decay entirely, the pulse propagating through the jet 
experiences defocusing due to electrons freed during the previous passes. These electrons have 
been liberated from their parent atom for a sufficiently long time, and have equilibrated into 



a true plasma (note that the free-electron states included in the quantum model that drives 
the pulse propagator are of different nature: they did not have enough time to thermalize, and 
their interactions is mainly with the nearby parent ion). Therefore, the influence of the plasma 
remnant can be included within the linear chromatic properties of the gas. Here we also assume 
that the diffusion resulted in a nearly homogeneous spatial distribution of free electrons, and 
the corresponding modification of the refractive index is constant over the cross-section of the 
beam. Of course, the spectral nature of our pulse solver allows to endow this susceptibility 
contribution with the correct frequency dependence, #(&)) ~ ~ 6} pi as / 6)2 ■ 

Figure [9] shows simulated and experimental spectra of harmonics from 7th to 15th. While 
details of relative strength of harmonics 9th to 15th depend on the exact placement of the beam 
focus with respect to the center of the jet (compare the two panels on the left), and on the 
absorption included in our simulation (compare full and dashed lines), we see that harmonics 
11 and 13 are the most pronounced, and that harmonic power starts to decrease at harmonics 
15th (there are many more harmonics generated, but are not discernible on the linear scale of 
this figure). This is compatible with the spectrum recorded in the experiment. On the other 
hand, the simulated 7th harmonic seems to be too strong. We think that this is mainly due to 
the fact that the model susceptibility of the gas has no absorption at this wavelength (this is 
due to limited frequency extent of the available data, see the red line in Fig. Q]). The effect of 
medium absorption is indeed clearly visible in the 9th harmonic where our model absorption 
data exhibit a maximum. These results make it thus evident that it is important to obtain as 
realistic as possible a data set for both the index of refraction and absorption of the gas. 
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Fig. 9. Harmonic spectrum calculated (left panels) for conditions reflecting those in the 
femtosecond enhancement cavity. Dashed and full line compare results for simulation with 
and without inclusion of losses. Left most panel: beam focus at the "entrance" into the gas 
jet. Middle panel: beam focus at the "exit" from the gas jet. Experimental spectrum shown 
in the right panel, with an account for the grating efficiency. 

We certainly do not expect that the simplistic one-dimensional quantum model put forward 
here could reproduce experimental results in a one-to-one manner. However, we see the above 
comparison as very encouraging — it shows that the model, especially when coupled with 
good-quality dispersion and absorption data for the gas, can serve as a practical tool to study 
the qualitative trends that govern the nonlinear interactions that span the frequency range from 
infrared to high harmonics. 



4. Conclusion 

In conclusion, we have presented a new computational approach for femtosecond pulse propa- 
gation in the transparency region of gases that permits full resolution in three space dimensions 
plus time while fully incorporating quantum coherent effects such as high-harmonic generation 
and strong-field ionization in a holistic fashion. The key innovation that makes this possible 
is the use of ID atomic model with delta-function potential to compute the nonlinear optical 
response due to ground state to continuum transitions, and the closed form solution for the non- 
linear response leads to significant reductions in computing time compared to a direct solution 
of the atomic Schrodinger equation. In spite of the approximation of a ID atomic model we con- 
tend that our computational approach will be of great utility as a research tool in a number of 
forefront areas. For example, in the field of femtosecond filamentation in gases an outstanding 
question is whether coherent light-matter interactions can play a significant role. The current 
models, based upon individual ingredients of the Kerr effect, multi-photon-ionization, and a 
Drude model for the plasma, are coming under increasing scrutiny. Our new approach answers 
to this need and provides a computationally viable means of performing simulations of fila- 
mentation which includes the quantum coherent atomic effects. Moreover, it eliminates the part 
of the filamentation model which is considered to be the weakest link, namely ionization and 
concomitant effects due to freed electrons. For this initial presentation we have chosen to illus- 
trate our computational approach using the example of HHG in a gas jet as this is a distinctly 
quantum coherent effect and also requites an enormous spectral width. In particular, using the 
example of near-threshold HHG in a Xenon gas jet we hope to have shown that our computa- 
tional approach yields results in keeping with general expectations, and that it is a useful tool 
to study qualitative dynamic trends in a completely self-consistent way. As a more concrete ex- 
ample we presented a qualitative comparison between our model and results from an in-house 
experiment on extreme ultraviolet generation in a femtosecond enhancement cavity. 



5. Appendix 

Here we summarize the formulas needed to evaluate the nonlinear component of the current 
induced in the system described by the time-dependent Schrodinger equation Q. 

An arbitrary time-dependent external field F(t) enters through the classical electron trajec- 
tory x c iit ) (here we assume that x is the direction of the optical field polarization) 



Pd{t) = -( F(r)dT , x d (t)= [ p c i(r)dr 
Jo Jo 



(6) 



Assuming that the system was initially in its ground state, an auxiliary quantity A (f ) is obtained 
first as a solution to the following integral equation: 



A{t) = \i/ R (-x c ,(t),t) + 



where the right-hand-side is 
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Having calculated A(t), the nonlinear (in the strength of the external field F) current is ex- 
pressed as a sum 

jW=jW+jM , (9) 



with the components listed below: 
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Details of the derivation and description of numerical implementation can be found in 
Ref. ll25l . An open-source implementation is also available upon request from the authors (M.K. 
or J.B.). 
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